# Replication file
rm(list=ls())
library(tidyverse)
library(dplyr)
library(writexl)
library(fixest)

############ Variable description 
# sp_bill: binary variable that takes 1 if a legislator (co)sponsors a bill
# demand_skilled:log demand for skilled foreign workers (standardized) (based on headquarters)
# wk_demand_skilled: log demand for skilled foreign workers using worksite addresses
# rep: takes 1 if a legistor is a Republican
# computer_occ: number of people in computer occupations (log)
# ba: number of people who graduate college (log standardized)
# unemp: number of people unemployed (log)
# median income: median income of a district (log)
# manufacture: number of people in manufacturing (log)
# agri: number of people in agriculture (log)
# citizen_naturalized:number of citizens naturalized (standardized)
# non_citizen:number of non-citizens (standardized)
# black:number of African Americans
# hipanic: number of hispanic
# asian:number of asian
# dem_share: two party democratic vote share
# number_firms: numebr of firms (headquartered) lobbying for immigration bills 
# party_code:100 democrat 200 republican
# nominate_dim1: ideology score

############# Figure 1
# Load data
d <- read_csv("bill_by_president.csv")

d$president <- factor(d$president, levels = c( "Bill Clinton",  "George W. Bush",  "Barack Obama","Donald Trump"))  

president_names <- c(
  `Bill Clinton` = "Clinton",
  `George W. Bush` = "Bush",
  `Barack Obama` = "Obama",
  `Donald Trump` = "Trump"
)

ggplot(d,                         
       aes(x = type,
           y = count,
           fill = bill_characteristics)) + 
  ylab("The number of bills introduced")+xlab("Bill type")+labs(fill="Characteristics")+
  scale_fill_hue(labels = c("Restricting", "Expanding", "Zero-sum")  )+
  geom_bar(stat = "identity",
           position = "stack") +
  facet_grid(~ president , labeller = as_labeller(president_names))+
  theme(legend.position = "top") + 
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank()) + 
  theme(legend.title = element_text( size=6.5), legend.text=element_text(size=6.5))

############### Figure 2
# Load data 
df <- read_csv("data_all.csv")

df116 <- df %>% 
  filter(congress_number==116) %>% 
  group_by(bill_name) %>% 
  mutate(dem_ratio = sum(sp_bill == 1 & party_code==100)/sum(party_code==100)) %>% 
  mutate(rep_ratio = sum(sp_bill == 1 & party_code==200)/sum(party_code==200)) %>% 
  mutate(distance = abs(dem_ratio-rep_ratio)) %>% 
  mutate(sd = sd(nominate_dim1, na.rm=TRUE))

unique116 <- df116[!duplicated(df116$distance), ]

#115th
df115 <- df %>% 
  filter(congress_number==115) %>% 
  group_by(bill_name) %>% 
  mutate(dem_ratio = sum(sp_bill == 1 & party_code==100)/sum(party_code==100)) %>% 
  mutate(rep_ratio = sum(sp_bill == 1 & party_code==200)/sum(party_code==200)) %>% 
  mutate(distance = abs(dem_ratio-rep_ratio)) %>% 
  mutate(sd = sd(nominate_dim1, na.rm=TRUE))

unique115 <- df115[!duplicated(df115$distance), ]

#114th
df114 <- df %>% 
  filter(congress_number==114) %>% 
  group_by(bill_name) %>% 
  mutate(dem_ratio = sum(sp_bill == 1 & party_code==100)/sum(party_code==100)) %>% 
  mutate(rep_ratio = sum(sp_bill == 1 & party_code==200)/sum(party_code==200)) %>% 
  mutate(distance = abs(dem_ratio-rep_ratio)) %>% 
  mutate(sd = sd(nominate_dim1, na.rm=TRUE))

unique114 <- df114[!duplicated(df114$distance), ]

#113th
df113 <- df %>%
  filter(congress_number==113) %>% 
  group_by(bill_name) %>% 
  mutate(dem_ratio = sum(sp_bill == 1 & party_code==100)/sum(party_code==100)) %>% 
  mutate(rep_ratio = sum(sp_bill == 1 & party_code==200)/sum(party_code==200)) %>% 
  mutate(distance = abs(dem_ratio-rep_ratio)) %>% 
  mutate(sd = sd(nominate_dim1, na.rm=TRUE))

unique113 <- df113[!duplicated(df113$distance), ]

#112
df112 <- df %>%
  filter(congress_number==112) %>% 
  group_by(bill_name) %>% 
  mutate(dem_ratio = sum(sp_bill == 1 & party_code==100)/sum(party_code==100)) %>% 
  mutate(rep_ratio = sum(sp_bill == 1 & party_code==200)/sum(party_code==200)) %>% 
  mutate(distance = abs(dem_ratio-rep_ratio)) %>% 
  mutate(sd = sd(nominate_dim1, na.rm=TRUE))

unique112 <- df112[!duplicated(df112$distance), ]

#111
df111 <- df %>%
  filter(congress_number==111) %>% 
  group_by(bill_name) %>% 
  mutate(dem_ratio = sum(sp_bill == 1 & party_code==100)/sum(party_code==100)) %>% 
  mutate(rep_ratio = sum(sp_bill == 1 & party_code==200)/sum(party_code==200)) %>% 
  mutate(distance = abs(dem_ratio-rep_ratio)) %>% 
  mutate(sd = sd(nominate_dim1, na.rm=TRUE))

unique111 <- df111[!duplicated(df111$distance), ]

#110
df110 <- df %>%
  filter(congress_number==110) %>% 
  group_by(bill_name) %>% 
  mutate(dem_ratio = sum(sp_bill == 1 & party_code==100)/sum(party_code==100)) %>% 
  mutate(rep_ratio = sum(sp_bill == 1 & party_code==200)/sum(party_code==200)) %>% 
  mutate(distance = abs(dem_ratio-rep_ratio)) %>% 
  mutate(sd = sd(nominate_dim1, na.rm=TRUE))

unique110 <- df110[!duplicated(df110$distance), ]

#109
df109 <- df %>%
  filter(congress_number==109) %>% 
  group_by(bill_name) %>% 
  mutate(dem_ratio = sum(sp_bill == 1 & party_code==100)/sum(party_code==100)) %>% 
  mutate(rep_ratio = sum(sp_bill == 1 & party_code==200)/sum(party_code==200)) %>% 
  mutate(distance = abs(dem_ratio-rep_ratio)) %>% 
  mutate(sd = sd(nominate_dim1, na.rm=TRUE))

unique109 <- df109[!duplicated(df109$distance), ]

#108
df108 <- df108 %>% 
  filter(congress_number==108) %>% 
  group_by(bill_name) %>% 
  mutate(dem_ratio = sum(sp_bill == 1 & party_code==100)/sum(party_code==100)) %>% 
  mutate(rep_ratio = sum(sp_bill == 1 & party_code==200)/sum(party_code==200)) %>% 
  mutate(distance = abs(dem_ratio-rep_ratio)) %>% 
  mutate(sd = sd(nominate_dim1, na.rm=TRUE))

unique108 <- df108[!duplicated(df108$distance), ]

library(plyr)
data <- rbind.fill(unique116, unique115, unique114, unique113, unique112, unique111, unique110, unique109, unique108)
detach("package:plyr")

library(ggjoy)
ggplot(data,aes(x = distance, y = as.factor(congress_number))) +
  geom_joy(scale = 2,alpha = .5,rel_min_height = 0.01) + theme_joy()+ylab('Congress') +xlab('Distance (Rice index)')


################# Table 2
df_open <- df %>% 
  filter(bill_character=="open")

expanding_st_1 <- feglm(sp_bill~ demand_skilled+rep|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(expanding_st_1)

expanding_st_2 <- feglm(sp_bill~ demand_skilled*rep|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(expanding_st_2)

expanding_st_3 <- feglm(sp_bill~ demand_skilled*rep+ computer_occ +ba + unemp + medianincome + agri + manufacture |state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(expanding_st_3)

expanding_st_4 <- feglm(sp_bill~ demand_skilled*rep+ computer_occ + ba+ unemp + medianincome + agri + manufacture + citizen_naturalized + non_citizen + black + hispanic + asian + dem_share |state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(expanding_st_4)

etable(expanding_st_1, expanding_st_2, expanding_st_3, expanding_st_4,  label = "tab:etable", tex = TRUE)

############### Figure 4a
beta.hat <- coef(expanding_st_4)
vcov1 <- vcov(expanding_st_4)
z0 <- c(rep(0, 500), rep(1, 500))
dy.dx <- beta.hat["demand_skilled"] + beta.hat["demand_skilled:rep"]*z0
se.dy.dx <- sqrt(vcov1["demand_skilled", "demand_skilled"] + (z0^2)*vcov1["demand_skilled:rep", "demand_skilled:rep"] + 2*z0*vcov1["demand_skilled", "demand_skilled:rep"])
upr <- dy.dx + 1.96*se.dy.dx
lwr <- dy.dx - 1.96*se.dy.dx
ggplot(data=NULL, aes(x=z0, y=dy.dx)) +
  labs(x="Partisanship", y="Marginal Effects of demand|Partisanship") +
  geom_point(aes(z0, dy.dx),size = 1) +
  geom_hline(yintercept=0, size = 1, linetype=3) +
  geom_errorbar(aes(ymin=lwr, ymax=upr), alpha=0.3, width=.3)+
  scale_x_discrete(expand = c(0,0.2) , limits = c(0,1), labels=c('Democrat', 'Republican'))+theme_minimal()
#ggsave(file="congress_bill/expanding_marginal.jpg")


################ Figure 4b
# democrats
tmpdat <- df_open[, c("demand_skilled","unemp", "computer_occ", "ba","medianincome", "agri", "manufacture","non_citizen",
                      "black", "hispanic", "citizen_naturalized", "asian", "dem_share" )]

tmpdat <- na.omit(tmpdat)
tmpdat$state_abbrev <- "CA"
tmpdat$time <- "2017-01-10"
tmpdat$rep <- 0

jvalues <- with(tmpdat, seq(from = min(demand_skilled), to = max(demand_skilled),length.out = 100))

pp <- lapply(jvalues, function(j) {
  tmpdat$demand_skilled <- j
  predict(expanding_st_4, newdata = tmpdat, type = "response")
})

plotdat <- t(sapply(pp, function(x) {
  c(M = mean(x), quantile(x, c(0.25, 0.75)))
}))

plotdat <- as.data.frame(cbind(plotdat, jvalues))

colnames(plotdat) <- c("PredictedProbability", "Lower", "Upper", "Demand")
head(plotdat)

plotdat$rep <- 0

# data for republicans
rtmpdat <- df_open[, c("demand_skilled","unemp", "computer_occ", "ba","medianincome", "agri", "manufacture","non_citizen",
                       "black", "hispanic", "citizen_naturalized", "asian", "dem_share" )]

rtmpdat <- na.omit(rtmpdat)
rtmpdat$state_abbrev <- "CA"
rtmpdat$time <- "2017-01-10"
rtmpdat$rep <- 1

rjvalues <- with(rtmpdat, seq(from = min(demand_skilled), to = max(demand_skilled),length.out = 100))

rpp <- lapply(rjvalues, function(j) {
  rtmpdat$demand_skilled <- j
  predict(expanding_st_4, newdata = rtmpdat, type = "response")
})

rplotdat <- t(sapply(rpp, function(x) {
  c(M = mean(x), quantile(x, c(0.25, 0.75)))
}))

rplotdat <- as.data.frame(cbind(rplotdat, rjvalues))

colnames(rplotdat) <- c("PredictedProbability", "Lower", "Upper", "Demand")
head(rplotdat)

rplotdat$rep <- 1 

# combine the two datasets
library(plyr)
d <- rbind.fill(rplotdat, plotdat)
detach("package:plyr")
d$rep <- as.factor(d$rep)

ggplot(d, aes(x = Demand, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = rep), alpha = 0.2) +
  geom_line(aes(colour = rep), size=1)+ 
  theme_minimal() +theme(legend.position = "top")+
  guides(colour=guide_legend(title="Partisanship")) +
  guides(fill=guide_legend(title="Partisanship")) +
  scale_color_manual(values=c('blue','red'), labels = c("Democrat", "Republican"))+
  scale_fill_manual(values=c("blue", "red"), labels = c("Democrat", "Republican"))+
  ylab("Predicted probability") + xlab("log Demand for skilled foreign workers") 
#ggsave(file="congress_bill/expanding_prob.jpg")

################ Table 3
df_zero <- df %>% 
  filter(bill_character=="zero-sum")  

zerosum_st_1 <- feglm(sp_bill~ demand_skilled+rep|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(zerosum_st_1)

zerosum_st_2 <- feglm(sp_bill~ demand_skilled*rep|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(zerosum_st_2)

zerosum_st_3 <- feglm(sp_bill~ demand_skilled*rep+ computer_occ+ba+unemp+ medianincome+agri+manufacture|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(zerosum_st_3)

zerosum_st_4 <- feglm(sp_bill~ demand_skilled*rep+ computer_occ+ba+unemp+ medianincome+agri+manufacture +citizen_naturalized+ non_citizen+ black + hispanic + asian+dem_share|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(zerosum_st_4)

etable(zerosum_st_1, zerosum_st_2, zerosum_st_3, zerosum_st_4,  label = "tab:etable", tex = TRUE)

############### Figure 5a
beta.hat <- coef(zerosum_st_4)
vcov1 <- vcov(zerosum_st_4)
z0 <- c(rep(0, 500), rep(1, 500))
dy.dx <- beta.hat["demand_skilled"] + beta.hat["demand_skilled:rep"]*z0
se.dy.dx <- sqrt(vcov1["demand_skilled", "demand_skilled"] + (z0^2)*vcov1["demand_skilled:rep", "demand_skilled:rep"] + 2*z0*vcov1["demand_skilled", "demand_skilled:rep"])
upr <- dy.dx + 1.96*se.dy.dx
lwr <- dy.dx - 1.96*se.dy.dx
ggplot(data=NULL, aes(x=z0, y=dy.dx)) +
  labs(x="Partisanship", y="Marginal Effects of demand|Partisanship") +
  geom_point(aes(z0, dy.dx),size = 1) +
  geom_hline(yintercept=0, size = 1, linetype=3) +
  geom_errorbar(aes(ymin=lwr, ymax=upr), alpha=0.3, width=.3)+
  scale_x_discrete(expand = c(0,0.2) , limits = c(0,1), labels=c('Democrat', 'Republican'))+theme_minimal()
#ggsave(file="congress_bill/figure 5a.jpg")

################ Figure 5b
# democrats
zdtmpdat <- df_zero[, c("demand_skilled","unemp", "computer_occ", "ba","medianincome", "agri","manufacture","non_citizen",
                        "black", "hispanic", "citizen_naturalized", "asian", "dem_share")]

zdtmpdat <- na.omit(zdtmpdat)
zdtmpdat$state_abbrev <- "CA"
zdtmpdat$time <- "2017-02-16"
zdtmpdat$rep <- 0

zdjvalues <- with(zdtmpdat, seq(from = min(demand_skilled), to = max(demand_skilled),length.out = 100))

zdpp <- lapply(zdjvalues, function(j) {
  zdtmpdat$demand_skilled <- j
  predict(zerosum_st_4, newdata = zdtmpdat, type = "response")
})


zdplotdat <- t(sapply(zdpp, function(x) {
  c(M = mean(x), quantile(x, c(0.25, 0.75)))
}))

zdplotdat <- as.data.frame(cbind(zdplotdat, zdjvalues))

colnames(zdplotdat) <- c("PredictedProbability", "Lower", "Upper", "Demand")
head(zdplotdat)

zdplotdat$rep <- 0

##data for republicans
zrtmpdat <- df_zero[, c("demand_skilled","unemp", "computer_occ", "ba","medianincome", "agri","manufacture","non_citizen",
                        "black", "hispanic", "citizen_naturalized", "asian", "dem_share")]

zrtmpdat <- na.omit(zrtmpdat)
zrtmpdat$state_abbrev <- "CA"
zrtmpdat$time <- "2017-02-16"
zrtmpdat$rep <- 1

zrjvalues <- with(zrtmpdat, seq(from = min(demand_skilled), to = max(demand_skilled),length.out = 100))

zrpp <- lapply(zrjvalues, function(j) {
  zrtmpdat$demand_skilled <- j
  predict(zerosum_st_4, newdata = zrtmpdat, type = "response")
})

zrplotdat <- t(sapply(zrpp, function(x) {
  c(M = mean(x), quantile(x, c(0.25, 0.75)))
}))

zrplotdat <- as.data.frame(cbind(zrplotdat, zrjvalues))

# better names and show the first few rows
colnames(zrplotdat) <- c("PredictedProbability", "Lower", "Upper", "Demand")
head(zrplotdat)

zrplotdat$rep <- 1 

# combine the two datasets
library(plyr)
z <- rbind.fill(zdplotdat, zrplotdat)
detach("package:plyr")
class(z$rep)
z$rep <- as.factor(z$rep)

ggplot(z, aes(x = Demand, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = rep), alpha = 0.2) +
  geom_line(aes(colour = rep), size=1)+ 
  theme_minimal() +theme(legend.position = "top")+
  guides(colour=guide_legend(title="Partisanship")) +
  guides(fill=guide_legend(title="Partisanship")) +
  scale_color_manual(values=c('blue','red'), labels = c("Democrat", "Republican"))+
  scale_fill_manual(values=c("blue", "red"), labels = c("Democrat", "Republican"))+
  ylab("Predicted probability") + xlab("log Demand for skilled foreign workers") 
#ggsave(file="congress_bill/figure 5b.jpg")

################# Table 4
df_restrict <- df %>% 
  filter(bill_character=="reduce_number"|bill_character=="strong_monitor")

restrict_st_1 <- feglm(sp_bill~ demand_skilled+rep|state_abbrev + time, data = df_restrict, family=binomial(link="logit"))
summary(restrict_st_1)

restrict_st_2 <- feglm(sp_bill~ demand_skilled*rep|state_abbrev + time, data = df_restrict, family=binomial(link="logit"))
summary(restrict_st_2)

restrict_st_3 <- feglm(sp_bill~ demand_skilled*rep+computer_occ+ ba + unemp + medianincome+agri+manufacture|state_abbrev + time, data = df_restrict, family=binomial(link="logit"))
summary(restrict_st_3)

restrict_st_4 <- feglm(sp_bill~ demand_skilled*rep+ computer_occ +ba+ unemp + medianincome+agri+manufacture+citizen_naturalized + non_citizen + black + hispanic + asian+ dem_share|state_abbrev + time, data = df_restrict, family=binomial(link="logit"))
summary(restrict_st_4)

etable(restrict_st_1, restrict_st_2, restrict_st_3, restrict_st_4,  label = "tab:etable", tex = TRUE)

####################### Figure 6a
beta.hat <- coef(restrict_st_4)
vcov1 <- vcov(restrict_st_4)
z0 <- c(rep(0, 500), rep(1, 500))
dy.dx <- beta.hat["demand_skilled"] + beta.hat["demand_skilled:rep"]*z0
se.dy.dx <- sqrt(vcov1["demand_skilled", "demand_skilled"] + (z0^2)*vcov1["demand_skilled:rep", "demand_skilled:rep"] + 2*z0*vcov1["demand_skilled", "demand_skilled:rep"])
upr <- dy.dx + 1.6*se.dy.dx
lwr <- dy.dx - 1.6*se.dy.dx
ggplot(data=NULL, aes(x=z0, y=dy.dx)) +
  labs(x="Partisanship", y="Marginal Effects of demand|Partisanship") +
  geom_point(aes(z0, dy.dx),size = 1) +
  geom_hline(yintercept=0, size = 1, linetype=3) +
  geom_errorbar(aes(ymin=lwr, ymax=upr), alpha=0.3, width=.3)+
  scale_x_discrete(expand = c(0,0.2) , limits = c(0,1), labels=c("Democrat", "Republican"))+theme_minimal()
#ggsave(file="congress_bill/figure 6a.jpg")

################## Figure 6b
# democrats
zdtmpdat <- df_restrict[, c("demand_skilled","unemp", "computer_occ", "ba","medianincome", "agri","manufacture","non_citizen",
                            "black", "hispanic", "citizen_naturalized", "asian", "dem_share")]

zdtmpdat <- na.omit(zdtmpdat)
zdtmpdat$state_abbrev <- "CA"
zdtmpdat$time <- "2017-03-23"
zdtmpdat$rep <- 0

zdjvalues <- with(zdtmpdat, seq(from = min(demand_skilled), to = max(demand_skilled),length.out = 100))

zdpp <- lapply(zdjvalues, function(j) {
  zdtmpdat$demand_skilled <- j
  predict(restrict_st_4, newdata = zdtmpdat, type = "response")
})


zdplotdat <- t(sapply(zdpp, function(x) {
  c(M = mean(x), quantile(x, c(0.25, 0.75)))
}))

zdplotdat <- as.data.frame(cbind(zdplotdat, zdjvalues))

colnames(zdplotdat) <- c("PredictedProbability", "Lower", "Upper", "Demand")
head(zdplotdat)

zdplotdat$rep <- 0

##data for republicans
zrtmpdat <- df_restrict[, c("demand_skilled","unemp", "computer_occ", "ba","medianincome", "agri", "manufacture","non_citizen",
                            "black", "hispanic", "citizen_naturalized", "asian", "dem_share")]

zrtmpdat <- na.omit(zrtmpdat)
zrtmpdat$state_abbrev <- "CA"
zrtmpdat$time <- "2017-03-23"
zrtmpdat$rep <- 1

zrjvalues <- with(zrtmpdat, seq(from = min(demand_skilled), to = max(demand_skilled),length.out = 100))

zrpp <- lapply(zrjvalues, function(j) {
  zrtmpdat$demand_skilled <- j
  predict(restrict_st_4, newdata = zrtmpdat, type = "response")
})

zrplotdat <- t(sapply(zrpp, function(x) {
  c(M = mean(x), quantile(x, c(0.25, 0.75)))
}))

zrplotdat <- as.data.frame(cbind(zrplotdat, zrjvalues))

colnames(zrplotdat) <- c("PredictedProbability", "Lower", "Upper", "Demand")
head(zrplotdat)

zrplotdat$rep <- 1 

# combine the two datasets
library(plyr)
z <- rbind.fill(zdplotdat, zrplotdat)
detach("package:plyr")
class(z$rep)
z$rep <- as.factor(z$rep)

ggplot(z, aes(x = Demand, y = PredictedProbability)) + 
  geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = rep), alpha = 0.2) +
  geom_line(aes(colour = rep), size=1)+ 
  theme_minimal() +theme(legend.position = "top")+
  guides(colour=guide_legend(title="Partisanship")) +
  guides(fill=guide_legend(title="Partisanship")) +
  scale_color_manual(values=c('blue','red'), labels = c("Democrat", "Republican"))+
  scale_fill_manual(values=c("blue", "red"), labels = c("Democrat", "Republican"))+
  ylab("Predicted probability") + xlab("log Demand for skilled foreign workers") 
#ggsave(file="congress_bill/figure 6b.jpg")

############ Table 5
hq_expanding_1 <- feglm(sp_bill~ number_firms*rep+ computer_occ+ba+unemp+ medianincome + agri+manufacture+citizen_naturalized+non_citizen+black+hispanic+asian+dem_share|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(hq_expanding_1)

hq_zerosum_1 <- feglm(sp_bill~ number_firms*rep+ computer_occ+ba+unemp+ medianincome + agri+manufacture + citizen_naturalized+non_citizen+black+hispanic+asian+dem_share|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(hq_zerosum_1)

hq_restrict_1 <- feglm(sp_bill~ number_firms*rep+ computer_occ+ba+unemp+ medianincome + agri+manufacture + citizen_naturalized+non_citizen+black+hispanic+asian+dem_share|state_abbrev + time, data = df_restrict, family=binomial(link="logit"))
summary(hq_restrict_1)

etable(hq_expanding_1, hq_zerosum_1, hq_restrict_1,  label = "tab:etable", tex = TRUE)


########### Appendix ##########################
# Table A1
library(vtable)
sum_df <- df %>% 
  dplyr::select(sp_bill, demand_skilled, rep, computer_occ, ba,unemp, medianincome, agri, manufacture,dem_share,citizen_naturalized, non_citizen, black, hispanic, asian)

st(sum_df, out="latex")

# Table A2
w_expanding_st_1 <- feglm(sp_bill~ wk_demand_skilled+rep|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(w_expanding_st_1)

w_expanding_st_2 <- feglm(sp_bill~ wk_demand_skilled*rep|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(w_expanding_st_2)

w_expanding_st_3 <- feglm(sp_bill~ wk_demand_skilled*rep+ computer_occ+ba+unemp  + medianincome+ agri + manufacture |state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(w_expanding_st_3)

w_expanding_st_4 <- feglm(sp_bill~ wk_demand_skilled*rep+ computer_occ+ba+unemp+ medianincome + agri+manufacture+ citizen_naturalized+non_citizen+black+hispanic+asian+dem_share|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(w_expanding_st_4)

etable(w_expanding_st_1, w_expanding_st_2, w_expanding_st_3, w_expanding_st_4,  label = "tab:etable", tex = TRUE)

# Table A3
w_zerosum_st_1 <- feglm(sp_bill~ wk_demand_skilled+rep|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(w_zerosum_st_1)

w_zerosum_st_2 <- feglm(sp_bill~ wk_demand_skilled*rep|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(w_zerosum_st_2)

w_zerosum_st_3 <- feglm(sp_bill~ wk_demand_skilled*rep+ computer_occ+ba+unemp  + medianincome+ agri + manufacture |state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(w_zerosum_st_3)

w_zerosum_st_4 <- feglm(sp_bill~ wk_demand_skilled*rep+ computer_occ+ba+unemp  + medianincome+ agri + manufacture + citizen_naturalized+non_citizen+black+hispanic+asian+dem_share|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(w_zerosum_st_4)

etable(w_zerosum_st_1, w_zerosum_st_2, w_zerosum_st_3, w_zerosum_st_4,  label = "tab:etable", tex = TRUE)

# Table A4
w_restrict_st_1 <- feglm(sp_bill~ wk_demand_skilled+rep|state_abbrev + time, data = df_restrict, family=binomial(link="logit"))
summary(w_restrict_st_1)

w_restrict_st_2 <- feglm(sp_bill~ wk_demand_skilled*rep|state_abbrev + time, data = df_restrict, family=binomial(link="logit"))
summary(w_restrict_st_2)

w_restrict_st_3 <- feglm(sp_bill~ wk_demand_skilled*rep+ computer_occ+ba+unemp  + medianincome+ agri + manufacture |state_abbrev + time, data = df_restrict, family=binomial(link="logit"))
summary(w_restrict_st_3)

w_restrict_st_4 <- feglm(sp_bill~ wk_demand_skilled*rep+ computer_occ+ba+unemp  + medianincome+ agri + manufacture + citizen_naturalized+non_citizen+black+hispanic+asian+dem_share|state_abbrev + time, data = df_restrict, family=binomial(link="logit"))
summary(w_restrict_st_4)

etable(w_restrict_st_1, w_restrict_st_2, w_restrict_st_3, w_restrict_st_4,  label = "tab:etable", tex = TRUE)

# Table A5 
ba1 <- feglm(sp_bill~ ba+rep|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(ba1)
ba2 <- feglm(sp_bill~ ba*rep|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(ba2)
ba3 <- feglm(sp_bill~ ba*rep+ unemp+ medianincome + agri+manufacture |state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(ba3)
ba4 <- feglm(sp_bill~ ba*rep+ unemp+ medianincome +agri+ manufacture +citizen_naturalized+non_citizen+black+hispanic+asian+dem_share|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(ba4)

compare_ba_5 <- feglm(sp_bill~ demand_skilled*rep+ba+ unemp+ medianincome +agri+manufacture+citizen_naturalized+non_citizen+black+hispanic+asian+dem_share|state_abbrev + time, data = df_open, family=binomial(link="logit"))
summary(compare_ba_5)

etable(ba1, ba2, ba3, ba4, compare_ba_5, label = "tab:etable", tex = TRUE)

# Table A6 
zerosum_ba_1 <- feglm(sp_bill~ demand_skilled+rep|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(zerosum_ba_1)

zerosum_ba_2 <- feglm(sp_bill~ demand_skilled*rep|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(zerosum_ba_2)

zerosum_ba_3 <- feglm(sp_bill~ demand_skilled*rep+ computer_occ+ba+unemp+ medianincome+agri+manufacture|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(zerosum_ba_3)

zerosum_ba_4 <- feglm(sp_bill~ demand_skilled*rep+ computer_occ+ba+unemp+ medianincome+agri+manufacture +citizen_naturalized+non_citizen+black+hispanic+asian+dem_share|state_abbrev + time, data = df_zero, family=binomial(link="logit"))
summary(zerosum_ba_4)

etable(zerosum_ba_1, zerosum_ba_2, zerosum_ba_3, zerosum_ba_4,  label = "tab:etable", tex = TRUE)

# Table A7 
robust_expanding_6 <- glm(sp_bill~ demand_skilled*rep+computer_occ+ba+unemp + medianincome +agri+ manufacture + dem_share +citizen_naturalized+non_citizen+black+hispanic+asian+factor(district) + factor(time)-1, data = df_open, family="binomial")
summary(robust_expanding_6)

robust_zerosum_6 <- glm(sp_bill~ demand_skilled*rep+computer_occ+ba+unemp + medianincome +agri+ manufacture + dem_share+citizen_naturalized+non_citizen+black+hispanic+asian+factor(district) + factor(time)-1, data = df_zero, family="binomial")
summary(robust_zerosum_6)

robust_restrict_6 <- glm(sp_bill~ demand_skilled*rep+computer_occ+ba+unemp + medianincome +agri+ manufacture + dem_share +citizen_naturalized+non_citizen+black+hispanic+asian+factor(district) + factor(time)-1, data = df_restrict, family="binomial")
summary(robust_restrict_6)

library(stargazer)
stargazer(robust_expanding_6, robust_zerosum_6, robust_restrict_6, type="latex", omit=c("district", "time"))

# Figure A3
df$rep <- as.factor(df$rep)
ggplot(df, aes(demand_skilled, fill = rep)) + 
  geom_histogram(alpha = 0.5, aes(y = after_stat(density)), position = 'identity') + theme_minimal() +
  theme(legend.position = "top") + labs(x="log Demand for skilled labor", y="Density",fill="Republican") 


# Figure A4
library(tidyverse)
library(broom)
library(dotwhisker)
t <- (tidy(expanding_st_4, exponentiate = FALSE, conf.int = TRUE)
      |> filter(stringr::str_detect(term, "^log|rep|computer|medianincome|agri|manufacture|dem"))
)

dwplot(t, vline = geom_vline(
  xintercept = 0,
  #colour = "grey60",
  linetype = 2),
  vars_order = c("demand_skilled", "rep", "demand_skilled:rep", "computer_occ", "ba","unemp", 
                 "medianincome",  "agri", "manufacture" ,"citizen_naturalized", "non_citizen",
                 "black", "hispanic",   "asian" , "dem_share")
)%>% 
  relabel_predictors(
    c(demand_skilled = "Demand",
      rep = "Republican",
      "demand_skilled:rep" = "Demand*Republican",
      computer_occ= "Computer occupations",
      ba="College graduates",
      unemp = "Unemployment",
      medianincome = "Median income",
      agri = "Agriculture", 
      manufacture = "Manufacturing",
      dem_share = "Dem vote share",
      citizen_naturalized = "Nuturalized citizens",
      non_citizen ="Non-citizens", 
      black = "African American",
      hispanic = "Hispanic", 
      asian = "Asian")) + 
  xlab("Coefficient estimate") + ylab("") +theme_minimal() + theme(legend.position="none")

# 
library(tidyverse)
library(broom)
library(dotwhisker)
z <- (tidy(zerosum_st_4, exponentiate = FALSE, conf.int = TRUE)
      |> filter(stringr::str_detect(term, "^log|rep|computer|medianincome|agri|manufacture|dem"))
)

dwplot(z, vline = geom_vline(
  xintercept = 0,
  #colour = "grey60",
  linetype = 2),
  vars_order = c("demand_skilled", "rep", "demand_skilled:rep", "computer_occ", "ba","unemp", 
                 "medianincome",  "agri", "manufacture"  ,"citizen_naturalized", "non_citizen",
                 "black", "hispanic",   "asian" , "dem_share" )
)%>% 
  relabel_predictors(
    c(demand_skilled = "Demand",
      rep = "Republican",
      "demand_skilled:rep" = "Demand*Republican",
      computer_occ= "Computer occupations",
      ba="College graduates",
      unemp = "Unemployment",
      medianincome = "Median income",
      agri = "Agriculture", 
      manufacture = "Manufacturing",
      dem_share = "Dem vote share",
      citizen_naturalized = "Nuturalized citizens",
      non_citizen ="Non-citizens", 
      black = "African American",
      hispanic = "Hispanic", 
      asian = "Asian")) +  
  xlab("Coefficient estimate") + ylab("") +theme_minimal() + theme(legend.position="none")


